This script generates plots comparing grassland bird population trends in Crow-Hassan park with regional trends from the BBS.
knitr::opts_chunk$set(echo = TRUE)
library(here)
library(tidyverse)
library(auk)
source(here("Scripts", "1_Data_Import.R"))
# source(here("Scripts", "2_Spatial_Data_Import.R"))
# source(here("Scripts", "3_Data_Cleaning.R"))
# source(here("Scripts", "4_Covariate_Data_Prep.R"))
source(here("Scripts", "5_Prep_Species.R"))# Which park estimates to use? Defined using the tag used to run the bird_trends_species script.
tag_chbbs <- "20231128_plus1_to_counts"
# Which regional estimates to use?
BBS_indices <- readRDS(here("Results", "993_BBS_fromCWS", "Indices", "00_all_sp_MN23.RDS"))
BBS_trends <- read_csv(here("Data", "Wildlife", "BBS", "BBS_CWS", "All_2021_BBS_trends.csv")) %>%
filter(Region == "US-MN-23")
# Create directory to save outputs
mainDir <- here("Results", "979_BBS_CHBB_comparison")
dir.create(file.path(mainDir, tag_chbbs), showWarnings = TRUE)
outputDir <- here("Results", "979_BBS_CHBB_comparison", tag_chbbs)
subdirectories <- c("Plots")
lapply(file.path(mainDir, tag_chbbs, subdirectories), function(x) if (!dir.exists(x)) dir.create(x))# Load park abundance estimates
chbbs_ests <- list.files(here("Results", "991_bird_trends_species", tag_chbbs, "annual_n_estimates"), pattern = "annual_N_ests.csv", full.names = TRUE)
chbbs_ests <- lapply(chbbs_ests, read_csv)
chbbs_ests_names <- list.files(here("Results", "991_bird_trends_species", tag_chbbs, "annual_n_estimates"), pattern = "annual_N_ests.csv", full.names = FALSE)
chbbs_ests_names <- gsub(pattern = "_annual_N_ests.csv", replacement = "", x = chbbs_ests_names)
names(chbbs_ests) <- chbbs_ests_names
chbbs_ests <- bind_rows(chbbs_ests, .id = "sp") %>%
select(sp, year, N_est_mean, N_est_LCI_95, N_est_UCI_95) %>%
rename(
Year = year,
Index = N_est_mean,
Index_q_0.025 = N_est_LCI_95,
Index_q_0.975 = N_est_UCI_95
) %>%
mutate(source = "chbbs")
# join eBird taxonomy and get list of park species
chbbs_ests <- chbbs_ests %>%
left_join(bird.tblBirdSpecies.ebird %>% select(BirdCode, common_name, scientific_name, taxon_order), by = c("sp" = "BirdCode")) %>%
mutate(species = common_name)
sp_chbbs <- unique(chbbs_ests$common_name)dat <- bind_rows(chbbs_ests, bbs_ests) %>%
mutate(source = factor(source, levels = c("chbbs", "bbs")))
adjust_chbbs_by_1 <- TRUE
if (adjust_chbbs_by_1) {
dat2 <- dat %>%
mutate(
Index_q_0.025 = ifelse(source == "chbbs", Index_q_0.025 - 1, Index_q_0.025),
Index_q_0.975 = ifelse(source == "chbbs", Index_q_0.975 - 1, Index_q_0.975),
Index = ifelse(source == "chbbs", Index - 1, Index)
) %>%
mutate(
Index_q_0.025 = ifelse(Index_q_0.025 < 0, 0, Index_q_0.025),
Index_q_0.975 = ifelse(Index_q_0.975 < 0, 0, Index_q_0.975),
Index = ifelse(Index < 0, 0, Index)
)
} else {
dat2 <- dat
}
sp_shared <- intersect(sp_bbs, sp_chbbs)
# Note species not found in one dataset or the other
# (sp_not_in_chbbs <- setdiff(sp_bbs, sp_chbbs)) # species in bbs but not chbbs
# (sp_not_in_bbs <- setdiff(sp_chbbs, sp_bbs)) # species in chbbs but not bbsgenerate_comp_plots <- function(focal_species) {
# Annual estimate plots
plot_chbbs_prop <- dat2 %>%
filter(species == focal_species) %>%
ggplot() +
ggthemes::theme_wsj(color = "grey") +
geom_ribbon(aes(x = Year, ymin = Index_q_0.025, ymax = Index_q_0.975, fill = source), alpha = 0.2) +
geom_line(aes(x = Year, y = Index, color = source, group = source)) +
geom_point(aes(x = Year, y = Index, color = source)) +
scale_color_manual(values = c("#24537D", "#5f9a90")) +
scale_fill_manual(values = c("#24537D", "#5f9a90")) +
facet_wrap(~source,
scales = "free_y", strip.position = "top", ncol = 1,
labeller = as_labeller(c(chbbs = "Total count (TRPD Prairie Transects)", bbs = "Mean count per route (MN-23 Breeding Bird Survey)"))
) +
ylab(NULL) +
theme(
strip.background = element_blank(),
strip.placement = "outside",
strip.text = element_text(hjust = 0, face = "italic", size = 10),
axis.text = element_text(face = "bold", size = 10),
plot.title = element_text(face = "bold", size = 12, family = NA),
legend.position = "none"
) +
labs()
# Percent change plot
# This calculates percent change in Index relative to the first shared year between the two data sets for each species
plot_percent_change <- dat %>%
filter(species == focal_species) %>%
group_by(species) %>%
# Find the first shared year between Methods for each species
mutate(First_Shared_Year = min(intersect(Year[source == "bbs"], Year[source == "chbbs"]))) %>%
ungroup() %>%
group_by(species, source) %>%
# Calculate percent change relative to the first shared year
mutate(Percent_Change = (Index - Index[Year == First_Shared_Year]) / Index[Year == First_Shared_Year] * 100) %>%
# Plot data
ggplot() +
ggthemes::theme_wsj(color = "grey") +
ggthemes::scale_color_wsj() +
geom_hline(yintercept = 0, color = "black", linewidth = 1.2, alpha = .5) +
geom_point(aes(x = Year, y = Percent_Change, color = source)) +
geom_line(aes(x = Year, y = Percent_Change, color = source, group = source)) +
facet_wrap(~species, scales = "free_y", strip.position = "top") +
scale_y_continuous(labels = function(x) paste0(x, "%")) +
scale_color_manual(values = c("#24537D", "#5f9a90")) +
scale_fill_manual(values = c("#24537D", "#5f9a90")) +
theme(
strip.background = element_blank(),
strip.placement = "outside",
# strip.text = element_text(hjust = 0, face = "italic", size = 10),
strip.text = element_blank(),
axis.text = element_text(face = "bold", size = 10),
plot.title = element_text(face = "bold", size = 12, family = NA),
plot.subtitle = element_text(hjust = 0, family = NA, face = "italic", size = 10),
legend.position = "none"
) +
labs( # title = focal_species,
subtitle = "Percent change relative to first shared year"
)
# Combine plots
plot_row <- cowplot::plot_grid(plot_chbbs_prop, plot_percent_change)
# Format and place title
title <- cowplot::ggdraw() +
cowplot::draw_label(
focal_species,
fontface = "bold",
x = 0,
hjust = 0
) +
theme(
# add margin on the left of the drawing canvas,
# so title is aligned with left edge of first plot
plot.margin = margin(0, 0, 0, 7)
)
# Layout
layout <- cowplot::plot_grid(
title, plot_row,
ncol = 1,
# rel_heights values control vertical title margins
rel_heights = c(0.05, 1)
)
# Save plot
ggsave(filename = paste0(focal_species, "_comp.png"), plot = layout, path = here("Results", "979_BBS_CHBB_comparison", tag_chbbs, "Plots"), width = 9, height = 5)
}
# Apply the function to all shared species
lapply(sp_shared, generate_comp_plots)files <- list.files(path = here("Results", "979_BBS_CHBB_comparison", tag_chbbs, "Plots"), pattern = "_comp.png", full.names = TRUE)
files_names <- list.files(path = here("Results", "979_BBS_CHBB_comparison", tag_chbbs, "Plots"), pattern = "_comp.png", full.names = FALSE)
files_names <- gsub("\\_comp.png", "", files_names)
for (i in 1:length(files)) {
cat(paste0("#### ", files_names[i], "\n"))
cat(paste0("\n\n"))
}# BBS trend table
BBS_trends$species[which(BBS_trends$species == "Red-tailed Hawk (all forms)")] <- "Red-tailed Hawk"
BBS_trends$species[which(BBS_trends$species == "Northern Flicker (all forms)")] <- "Northern Flicker"
BBS_trends_summary <- BBS_trends %>%
filter(Trend_Time == "Long-term") %>%
mutate(source = "bbs") %>%
mutate(
Trend_summary = ifelse(Trend_Q0.025 < 0 & Trend_Q0.975 < 0, "Decreasing",
ifelse(Trend_Q0.025 <= 0 & Trend_Q0.975 >= 0, "Stable or uncertain", "Increasing")
),
converged = 1,
converged = as.factor(converged)
) %>%
left_join(bird.tblBirdSpecies.ebird %>% select(BirdCode, common_name, taxon_order), by = c("species" = "common_name")) %>%
rename(common_name = species) %>%
select(common_name, Start_year, End_year, Relative_Abundance, source, Trend, Trend_Q0.025, Trend_Q0.975, Trend_summary, converged, Percent_Change, Percent_Change_Q0.025, Percent_Change_Q0.975)
# CHBBS trend table
# Get max rhat of model for each species
species_level_estimates <- list.files(here("Results", "991_bird_trends_species", tag_chbbs, "species_level_estimates"), pattern = "species_level_estimates.csv", full.names = TRUE) %>%
lapply(read.csv) %>%
bind_rows() %>%
select(-X)
rhat_summary <- data.frame(
sp = species_level_estimates$sp,
max_rhat = species_level_estimates$max.rhat.run
)
# Summarize effect of year for each species (temporal trend)
CHBBS_trends_summary <- species_level_estimates %>%
select(sp, r0_mean, r0_sd, r0_LCI_95, r0_UCI_95, f_r0) %>%
left_join(rhat_summary, by = "sp")
CHBBS_trends_summary <- CHBBS_trends_summary %>%
mutate(
Trend_summary = ifelse(r0_LCI_95 < 0 & r0_UCI_95 < 0, "Decreasing",
ifelse(r0_LCI_95 <= 0 & r0_UCI_95 >= 0, "Stable or uncertain", "Increasing")
),
converged = ifelse(max_rhat >= 1.1, 0, 1),
converged = as.factor(converged)
) %>%
left_join(bird.tblBirdSpecies.ebird %>% select(BirdCode, common_name, taxon_order), by = c("sp" = "BirdCode"))
CHBBS_trends_summary$common_name <- reorder(CHBBS_trends_summary$common_name, CHBBS_trends_summary$r0_mean)
CHBBS_trends_summary$Trend_summary <- factor(CHBBS_trends_summary$Trend_summary, levels = c("Decreasing", "Stable or uncertain", "Increasing"))
CHBBS_trends_summary <- CHBBS_trends_summary %>%
mutate(
BirdCode = sp,
source = "chbbs"
) %>%
rename(
"Trend" = r0_mean,
"Trend_Q0.025" = r0_LCI_95,
"Trend_Q0.975" = r0_UCI_95
) %>%
select(common_name, Trend, Trend_Q0.025, Trend_Q0.975, source, Trend_summary, converged, BirdCode, taxon_order)
# Combined trends tables for plotting
dat_trends <- bind_rows(CHBBS_trends_summary, BBS_trends_summary)
dat_trends$common_name[which(dat_trends$common_name == "Common Snipe")] <- "Wilson's Snipe"
dat_trends_pivot <- dat_trends %>%
pivot_wider(id_cols = common_name, names_from = source, values_from = c(starts_with("Trend"), converged)) %>%
left_join(bird.tblBirdSpecies.ebird, by = "common_name") %>%
mutate(Trend_bbs = Trend_bbs / 100)
sp_detected_above_threshold <- readRDS(here("Results", "991_bird_trends_species", tag_chbbs, "sp_detected_above_threshold.RDS"))
species_to_run_above_threshold <- readRDS(here("Results", "991_bird_trends_species", tag_chbbs, "species_to_run_above_threshold.RDS"))# Generate scatter plots comparing regional and park trends
# Remove species observed <5 times before plotting
scatter <- ggplot(
data = dat_trends_pivot %>%
filter(BirdCode %in% sp_detected_above_threshold),
aes(x = Trend_bbs, y = Trend_chbbs, label = BirdCode)
) +
geom_hline(yintercept = 0, color = "grey50") +
geom_vline(xintercept = 0, color = "grey50") +
# stat_smooth(method = "lm",
# se = FALSE,
# color = "grey",
# linewidth = .5) +
ggrepel::geom_label_repel(size = 2.4, box.padding = 0.02, label.padding = .075, arrow = arrow(length = unit(0.01, "npc"), type = "open"), max.time = 20, max.iter = 100000, force_pull = 5, seed = 1, max.overlaps = 20, fill = "grey95") +
# scale_x_continuous(limits = ggpmisc::symmetric_limits) +
# scale_y_continuous(limits = ggpmisc::symmetric_limits) +
ggh4x::coord_axes_inside(labels_inside = TRUE) +
# scale_x_continuous(limits = ~ c(-1, 1) * max(abs(.x))) +
# scale_y_continuous(limits = ~ c(-1, 1) * max(abs(.x))) +
scale_y_continuous(
labels = scales::percent,
position = "right"
) +
scale_x_continuous(
labels = scales::percent,
position = "top"
) +
labs(
x = "Regional trend (%/yr)",
y = "Park trend (%/yr)"
) +
theme_classic()
ggsave(plot = scatter, filename = "comp_scatter.png", path = here("Results", "979_BBS_CHBB_comparison", tag_chbbs, "Plots"), width = 8, height = 5)
scatterscatter_form <- ggplot(data = dat_trends_pivot %>% filter(BirdCode %in% sp_detected_above_threshold), aes(x = Trend_bbs, y = Trend_chbbs, label = BirdCode)) +
# geom_point(alpha = .8, color = "#32537d") +
geom_hline(yintercept = 0, color = "grey50") +
geom_vline(xintercept = 0, color = "grey50") +
# stat_smooth(method = "lm",
# se = FALSE,
# color = "grey",
# linewidth = .5) +
ggrepel::geom_label_repel(size = 2.4, box.padding = 0.02, label.padding = .075, arrow = arrow(length = unit(0.01, "npc"), type = "open"), max.time = 20, max.iter = 100000, force_pull = 5, seed = 1, max.overlaps = 20, fill = "grey95") +
# scale_x_continuous(limits = ggpmisc::symmetric_limits) +
# scale_y_continuous(limits = ggpmisc::symmetric_limits) +
scale_y_continuous(
limits = c(-.08, .08),
breaks = seq(-.08, .08, .02),
labels = scales::percent,
position = "right"
) +
scale_x_continuous(
limits = c(-.11, .11),
breaks = seq(-.12, .12, .02),
labels = scales::percent,
position = "top"
) +
ggh4x::coord_axes_inside(labels_inside = TRUE) +
# scale_x_continuous(limits = ~ c(-1, 1) * max(abs(.x))) +
# scale_y_continuous(limits = ~ c(-1, 1) * max(abs(.x))) +
labs(
x = "Regional trend (%/yr)",
y = "Park trend (%/yr)"
) +
theme_classic() #+
# theme(axis.text = element_blank())
ggsave(plot = scatter_form, filename = "comp_scatter_form.png", path = here("Results", "979_BBS_CHBB_comparison", tag_chbbs, "Plots"), width = 8, height = 5)
scatter_formscatter_form_fit <- ggplot(data = dat_trends_pivot %>% filter(BirdCode %in% sp_detected_above_threshold), aes(x = Trend_bbs, y = Trend_chbbs, label = BirdCode)) +
# geom_point(alpha = .8, color = "#32537d") +
geom_hline(yintercept = 0, color = "grey50") +
geom_vline(xintercept = 0, color = "grey50") +
stat_smooth(
method = "lm",
formula = "y ~ x",
se = FALSE,
color = "grey",
linewidth = .5
) +
ggrepel::geom_label_repel(size = 2.4, box.padding = 0.02, label.padding = .075, arrow = arrow(length = unit(0.01, "npc"), type = "open"), max.time = 20, max.iter = 100000, force_pull = 5, seed = 1, max.overlaps = 20, fill = "grey95") +
# scale_x_continuous(limits = ggpmisc::symmetric_limits) +
# scale_y_continuous(limits = ggpmisc::symmetric_limits) +
scale_y_continuous(
limits = c(-.08, .08),
breaks = seq(-.08, .08, .02),
labels = scales::percent,
position = "right"
) +
scale_x_continuous(
limits = c(-.11, .11),
breaks = seq(-.12, .12, .02),
labels = scales::percent,
position = "top"
) +
ggh4x::coord_axes_inside(labels_inside = TRUE) +
# scale_x_continuous(limits = ~ c(-1, 1) * max(abs(.x))) +
# scale_y_continuous(limits = ~ c(-1, 1) * max(abs(.x))) +
labs(
x = "Regional trend (%/yr)",
y = "Park trend (%/yr)"
) +
theme_classic() #+
# theme(axis.text = element_blank())
ggsave(plot = scatter_form_fit, filename = "comp_scatter_form_fit.png", path = here("Results", "979_BBS_CHBB_comparison", tag_chbbs, "Plots"), width = 8, height = 5)
scatter_form_fitView regression fit:
##
## Call:
## lm(formula = Trend_chbbs ~ Trend_bbs, data = dat_trends_pivot)
##
## Coefficients:
## (Intercept) Trend_bbs
## 0.007365 0.018995
View residuals (the larger the value the more the park is overperforming the region):
resids <- model$residuals
comp_test <- dat_trends_pivot[-c(model$na.action), ]
names(resids) <- comp_test$common_name
resids <- enframe(resids) %>%
arrange(desc(value))
residsGoal here is to get all potential prairie-associated species that breed in MN
habs <- read.csv(here("Data", "Wildlife", "Bird_habitat_associations_USFS.csv")) %>%
mutate(Common.name = str_replace(Common.name, "<92>", "'"))
# Get species with at least probable breeding in MN from last BBA
mnbba_sp <- read_excel(here("Data", "Wildlife", "BBA", "Species_Summary_Table.xlsx"), skip = 1, .name_repair = "universal")
mnbba_sp$Species[which(mnbba_sp$Species == "Le Conte's Sparrow")] <- "LeConte's Sparrow"
mnbba_sp$Species[which(mnbba_sp$Species == "Gray Jay")] <- "Canada Jay"
mn_prob_and_conf_breeding <- mnbba_sp$Species[which(mnbba_sp$Probable...4 > 0 & mnbba_sp$Confirmed...5 > 0)]
# This is the species ACAD or MN BBA list as prairie (plus Tree Swallow)
acad_prairie <- acad %>%
filter(`Primary Breeding Habitat` %in% c("Grasslands: Temperate", "Wetlands: Seasonally Wet Grasslands", "Open Country: Habitat Mosaic", "Open Country Aerial: Habitat Mosaic") |
`Secondary Breeding Habitat` %in% c("Grasslands: Temperate", "Wetlands: Seasonally Wet Grasslands", "Open Country: Habitat Mosaic", "Open Country Aerial: Habitat Mosaic") |
`Common Name` %in% mn_bba_habs_prairie$common_name |
`Common Name` %in% "Tree Swallow")
# This is the MN breeding species (ACAD version)
acad_mn <- acad %>%
filter(`Common Name` %in% mn_prob_and_conf_breeding)
# These are the species that are both prairie and MN breeders
acad_mn_prairie <- intersect(acad_prairie, acad_mn)
acad_mn_prairie$`Common Name`## [1] "Northern Shoveler" "Northern Pintail"
## [3] "Canvasback" "Redhead"
## [5] "Killdeer" "Ruddy Duck"
## [7] "Mourning Dove" "Gray Partridge"
## [9] "Ring-necked Pheasant" "Wilson's Phalarope"
## [11] "Upland Sandpiper" "Sharp-tailed Grouse"
## [13] "Greater Prairie-Chicken" "Chimney Swift"
## [15] "Common Nighthawk" "Sandhill Crane"
## [17] "Marbled Godwit" "Great Horned Owl"
## [19] "Wilson's Snipe" "Peregrine Falcon"
## [21] "Northern Harrier" "Turkey Vulture"
## [23] "Red-headed Woodpecker" "Swainson's Hawk"
## [25] "Red-tailed Hawk" "Short-eared Owl"
## [27] "American Kestrel" "Merlin"
## [29] "Willow Flycatcher" "Barn Swallow"
## [31] "Eastern Phoebe" "Western Kingbird"
## [33] "Eastern Kingbird" "Tree Swallow"
## [35] "Loggerhead Shrike" "Eastern Bluebird"
## [37] "Common Raven" "Horned Lark"
## [39] "Purple Martin" "Northern Rough-winged Swallow"
## [41] "Bank Swallow" "Cliff Swallow"
## [43] "Black-billed Magpie" "American Crow"
## [45] "Sedge Wren" "Common Yellowthroat"
## [47] "House Finch" "American Goldfinch"
## [49] "Brown-headed Cowbird" "Eastern Meadowlark"
## [51] "Western Meadowlark" "Brewer's Blackbird"
## [53] "Common Grackle" "Clay-colored Sparrow"
## [55] "Field Sparrow" "Vesper Sparrow"
## [57] "Lark Sparrow" "Savannah Sparrow"
## [59] "Grasshopper Sparrow" "Henslow's Sparrow"
## [61] "LeConte's Sparrow" "Song Sparrow"
## [63] "Northern Cardinal" "Blue Grosbeak"
## [65] "Dickcissel" "Bobolink"
## [67] "Red-winged Blackbird" "Gadwall"
## [69] "American Wigeon" "Blue-winged Teal"
# Trends for all prairie species in MN breeding bird atlas
compare <- data.frame(common_name = acad_mn_prairie$`Common Name`) %>%
left_join(dat_trends_pivot, by = "common_name") %>%
mutate(
Trend_chbbs = Trend_chbbs,
diff = Trend_chbbs - Trend_bbs
) %>%
select(common_name, BirdCode, diff, starts_with("Trend"), taxon_order)
# Trends for prairie species in Three Rivers species table (more inclusive)
compare2 <- data.frame(common_name = prairie_sp_table$`common_name`) %>%
left_join(dat_trends_pivot, by = "common_name") %>%
mutate(
Trend_chbbs = Trend_chbbs,
diff = Trend_chbbs - Trend_bbs
) %>%
select(common_name, BirdCode, diff, starts_with("Trend"), taxon_order)
# What prairie species are on Three Rivers table that aren't potential breeders according to BBA?
setdiff(compare2$common_name, compare$common_name)## [1] "Baird's Sparrow" "Chestnut-collared Longspur" "Nelson's Sparrow"
## [4] "Northern Bobwhite" "Northern Mockingbird" "Yellow Rail"
# Using the more-inclusive list for now
data <- compare2 %>%
arrange(desc(Trend_bbs)) %>%
mutate(common_name = factor(common_name, common_name))
# Plot shared species (those in both park and region)
lolliplot_shared_sp <- data %>%
filter(!(is.na(Trend_bbs) | is.na(Trend_chbbs))) %>%
ggplot() +
geom_hline(yintercept = 0) +
geom_segment(aes(x = common_name, xend = common_name, y = Trend_chbbs, yend = Trend_bbs), color = "grey") +
geom_point(aes(x = common_name, y = Trend_bbs), color = "grey", size = 3, alpha = .6) +
geom_point(aes(x = common_name, y = Trend_chbbs), color = "#24537D", size = 3, alpha = .6) +
coord_flip() +
# theme_ipsum() +
theme(
legend.position = "none",
) +
xlab("") +
ylab("Trend (% change per year)") +
scale_y_continuous(labels = scales::percent) +
ggthemes::theme_clean()
ggsave(lolliplot_shared_sp, filename = "lolliplot_shared_sp.png", path = here("Results", "979_BBS_CHBB_comparison", tag_chbbs, "Plots"), width = 8, height = 7)
lolliplot_shared_sp# Plot all BBS species, whether they have occurred in park surveys or not
lolliplot_all_bbs_sp <- data %>%
filter(!(is.na(Trend_bbs) & is.na(Trend_chbbs))) %>%
ggplot() +
geom_hline(yintercept = 0) +
geom_segment(aes(x = common_name, xend = common_name, y = Trend_chbbs, yend = Trend_bbs), color = "grey") +
geom_point(aes(x = common_name, y = Trend_bbs), color = "grey", size = 3, alpha = .6) +
geom_point(aes(x = common_name, y = Trend_chbbs), color = "#24537D", size = 3, alpha = .6) +
coord_flip() +
# theme_ipsum() +
theme(
legend.position = "none",
) +
xlab("") +
ylab("Trend (% change per year)") +
scale_y_continuous(labels = scales::percent) +
ggthemes::theme_clean()
ggsave(lolliplot_all_bbs_sp, filename = "lolliplot_all_bbs_sp.png", path = here("Results", "979_BBS_CHBB_comparison", tag_chbbs, "Plots"), width = 8, height = 7)
lolliplot_all_bbs_sp# Plot shared species- excluding those with just a few CH obvs
lolliplot_shared_sp_excluding_sparse_CH_sp <- data %>%
filter(BirdCode %in% species_to_run_above_threshold) %>%
filter(!(is.na(Trend_bbs) | is.na(Trend_chbbs))) %>%
ggplot() +
geom_hline(yintercept = 0) +
geom_segment(aes(x = common_name, xend = common_name, y = Trend_chbbs, yend = Trend_bbs), color = "grey") +
geom_point(aes(x = common_name, y = Trend_bbs), color = "grey", size = 3, alpha = .6) +
geom_point(aes(x = common_name, y = Trend_chbbs), color = "#24537D", size = 3, alpha = .6) +
coord_flip() +
# theme_ipsum() +
theme(
legend.position = "none",
) +
xlab("") +
ylab("Trend (% change per year)") +
scale_y_continuous(labels = scales::percent) +
ggthemes::theme_clean()
ggsave(lolliplot_shared_sp_excluding_sparse_CH_sp, filename = "lolliplot_shared_sp_excluding_sparse_CH_sp.png", path = here("Results", "979_BBS_CHBB_comparison", tag_chbbs, "Plots"), width = 8, height = 7)
lolliplot_shared_sp_excluding_sparse_CH_sp# Plot all BBS species, but species with just a few CH obs don't get a CH dot
lolliplot_all_bbs_sp_no_dot_for_sparse_CH_sp <- data %>%
filter(!(is.na(Trend_bbs) & is.na(Trend_chbbs))) %>%
mutate(Trend_chbbs = ifelse(BirdCode %in% species_to_run_above_threshold, Trend_chbbs, NA)) %>%
ggplot() +
geom_hline(yintercept = 0) +
geom_segment(aes(x = common_name, xend = common_name, y = Trend_chbbs, yend = Trend_bbs), color = "grey") +
geom_point(aes(x = common_name, y = Trend_bbs), color = "grey", size = 3, alpha = .6) +
geom_point(aes(x = common_name, y = Trend_chbbs), color = "#24537D", size = 3, alpha = .6) +
coord_flip() +
# theme_ipsum() +
theme(
legend.position = "none",
) +
xlab("") +
ylab("Trend (% change per year)") +
scale_y_continuous(labels = scales::percent) +
ggthemes::theme_clean()
ggsave(lolliplot_all_bbs_sp_no_dot_for_sparse_CH_sp, filename = "lolliplot_all_bbs_sp_no_dot_for_sparse_CH_sp.png", path = here("Results", "979_BBS_CHBB_comparison", tag_chbbs, "Plots"), width = 8, height = 7)
lolliplot_all_bbs_sp_no_dot_for_sparse_CH_spBy Sam Safran